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Abstract 

Motivated by the problem of solving the Einstein equations, we discuss 
high order finite difference discretizations of first order in time, second or- 
der in space hyperbolic systems. Particular attention is paid to the case 
when first order derivatives that can be identified with advection terms 
are approximated with non-centered finite difference operators. We first 
derive general properties of these discrete operators, then we extend a 
known result on numerical stability for such systems to general order of 
accuracy. As an application we analyze the shifted wave equation, in- 
cluding the behavior of the numerical phase and group speeds at different 
orders of approximations. Special attention is paid to when the use of 
off-centered schemes improves the accuracy over the centered schemes. 



1 Introduction 

Numerical discretization of first order hyperbolic systems of partial differential 
equations (PDEs) is greatly simplified by a result for the linear constant coeffi- 
cient case [1] : If the Cauchy problem is well-posed, then the semi-discrete prob- 
lem (only discretizing space and leaving time continuous) is stable when spatial 
derivatives are discretized with a centered finite difference operator (CFDO). 
Furthermore, when using simple Runge-Kutta methods [5] for time integration, 
for sufficiently small time step the fully discrete problem is also stable. 

Such a result does not hold in general for second order systems where first and 
second spatial derivatives appear [3J! In order to obtain a stable semi-discrete 
scheme, the second order system needs to have additional properties. In [3J, 
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which in the following we refer to as CHH, sufficient conditions for stability 
of the fully discrete problem were presented for such systems. Although these 
conditions arc valid for general order centered discretizations, this point has not 
been explicitly made. One of the results of this article is to make this statement 
clear, by closing a technical gap related to the boundedness of the lower order 
terms. The main focus of our work is on a detailed analysis of the numerical 
properties of discretizations where some hrst order derivatives are approximated 
with off-centcrcd finite difference operators and artificial dissipation is added to 
the equations. The motivation for choosing this more general situation comes 
from numerical relativity, where it is a common practice to off-center by one 
point the derivatives corresponding to the Lie advection terms. In numerical 
simulations of black holes using the BSSN formulation of the Einstein equations 
HH6], this procedure of off-centering was found to be essential even for sixth 
order schemes [7]. Numerical solutions of the Einstein equations are currently 
quickly expanding our knowledge about the astrophysics of compact binaries 
(see [5)15] for overviews on what has been achieved since the major breakthroughs 
in 2005 PUHUDi but a systematic understanding of the underlying numerical 
techniques has not yet been achieved. 

The shifted scalar wave equation serves as a simple but powerful model in 
numerical relativity [T5lU7j . The particular case with zero shift and flat back- 
ground (standard wave equation) has been extensively studied in [TS]- [H] and 
high order discretization methods have been proposed. In Sec. Owe introduce 
the shifted scalar wave equation as a first order in time, second order in space 
system, together with a summary of the well-posedness theory for mixed order 
systems. We show stability for our semi-discrete problem, independent of shift 
or dissipation terms, while the fully-discrete problem requires artificial dissipa- 
tion if more than one point is off-centered. Restricting to flat space in one space 
dimension, Courant limits and numerical phase and group speeds are computed 
and analyzed in detail. It is shown that increasing the off-centering reduces the 
Courant limit. However, by increasing the order of approximation while keep- 
ing the off-centering fixed, does not necessarily generate lower Courant limits. 
Regarding the numerical speeds, it is shown that indeed there are cases when 
off-centering improves the accuracy over the centered scheme. This fact is illus- 
trated also experimentally by the results of some simple numerical tests at the 
end of Section [5] 

Our analysis of the wave equation relies on certain properties of finite differ- 
ence operators, in particular on their behavior in Fourier space. We introduce 
these operators in Section|3]togcther with highlighting some relevant properties. 
Then, in Section [4] we address the stability method and follow in Section [5] with 
the analysis of the wave equation. Our results are summarized in Section [SJ 
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2 The shifted wave equation and first order in 
time second order in space hyperbolic systems 

The scalar wave equation in a d + 1-dimensional spacetime equipped with a 
Lorentzian metric g a p reads 

g ap d a dp0 = 0. (1) 

We assume a uniform time slicing for simplicity, g 00 = — 1 , and perform arf+l 
split introducing a positive definite d-metric ^ l = gi l + f3 l f3 3 , with i,j = l,d 
and a shift vector /3 l = g° l (see e.g. [2D])- The wave equation ([T]) then becomes 

d tt <t> = 2p%d t <t> + (y 1 - pi?) didj<t> . 

The mixed time-space derivatives lead to non-standard behavior as compared 
to the fiat space wave equation with zero shift, and much of the material below 
will be devoted to their treatment. 

We reduce the wave equation to a first order in time, second order in space 
form by introducing the variable K, in analogy with the York-ADM-system [2"T] 
(and other common representations of the Einstein equations), 

K = d t <t>- P'djQ 

which transforms the wave equation into the first order in time, second order in 
space system in the way most common in numerical relativity: 

d t <t> = ftdjQ + K, d t K = y jl d j i<t> + l3 j d j K. (2) 

Well-poscdness for the Cauchy problem for the system (JTJ is a standard textbook 
result both in the original second order form and for reduction to first order 
symmetric hyperbolic form. In the latter form standard theorems for numerical 
stability apply [1] . Here we investigate the numerical stability for the first order 
in time, second order in space system (J5|), using the methods presented in [3J. 
In this respect, the appropriate generalization of the shifted wave equation is a 
linear system of PDEs with constant coefficients of the form [3J : 

± V (t,x) = Pv(t,x), v=(U,V) T , 

with x G M. d , U : RxR d -> W, V : MxM d -> R 9 and 

/ A'dj + B C 

Note that the state vector v is split into two parts, U are those variables for 
which only first spatial derivatives appear, while second spatial derivatives of 
the ^-variables do enter the PDE. The well-posedness of the Cauchy problem 
for first order in time, second order in space systems of PDEs systems has been 




3 



clarified by |22fl2l)] . We will here recall the presentation in CHH, where the 
well-posedness of such systems of PDEs is discussed in close analogy with the 
issue of numerical stability. It is natural to consider 27r-periodic solutions and 
turn the analysis in Fourier space. 

In Fourier space the evolution problem reduces to a system of ordinary dif- 
ferential equations (ODEs) for the Fourier coefficients. By performing a first 
order reduction in Fourier space, it can be shown that well-posedness is not 
influenced by lower differential order terms, which we can therefore drop and 
consider the second order principal symbol constructed as 

/ iu:oA n C \ 
\ -oj$D nn iu} G n J ' 

where ojq = M, u> = (oji, . . . , 0J4) € M n = Ad^rij and rij = wjuIq 1 . 

It is shown in [3] that if there exists a matrix H(o>) = hr (cj) such that 
HP' + P'^H = and a positive constant K, such that K^ 1 !^ < H < KI^ 
(where I Uo = diag^/p, I q ]), then the problem is well-posed in the norms 

d 

El^U| 2 + |V| 2 , ||v|| 2 H = £ v f Hv- 

3=1 U)£Z d 

In [3] an analysis of numerical stability was performed in analogy with the 
proof of well-posedness as we have sketched it, and which we will extend to 
arbitrary approximation order in Section But, before going into stability 
analysis, we need to discuss some general properties of finite difference operators. 




3 Finite Difference Operators 

3.1 Construction and properties in one dimension 

Consider a mesh of equidistant points x v — vh, with veZ and h representing 
the grid spacing. Corresponding to the continuum vector function v : K — > C x 
■ ■ • x C we associate the grid vector function v by v : {x„, v = 0, ±1, ±2, ...}—> 
C x • • • x C and v v = v{x v ) = v(x„). 

Using 2n + 1 consecutive points, we want to construct the finite difference 
operator corresponding to the m-derivative. Let s £ {0, 1, . . . , n} be the offset 
of these points from symmetry with respect to the center, (s = for CFDO) and 
e the direction of off-centering (e = 1 for off-centering to the right, e = — 1 for 
off-centering to the left)!]] Then the finite difference operator to be constructed 
will be denoted £)( m ' n > s ' e ). It is a linear combination of shift operators of the 
form: 

n+es 

D ( m ,n,s,e) = h -m f mmk S k , (4) 

k=— n-\-es 

1 Though one can simplify the notation by dropping e and considering s £ {— n, . . . , n}, it 
will later turn out useful to separate the sign of s and its absolute value. 
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where S k be the shift operator by k points, S k v v = tv+fc- The weights f my n,s,e,k 
can be expressed as the coefficients of y k in the Taylor expansion of the function 



f m,n,s,e {y)=y n-e S 

around the point ?/o = 1 up to the order (y—yo) 2n (see appendix [5] for the proof). 
Using this procedure, one can deduce explicit expressions for the finite difference 
operators corresponding to the first and second derivative (the relations (|60p 
from appendix [A| . 

These expressions are fairly complicated, but they can be written in a more 
convenient form if we make use of the elementary finite difference operators: 

D±v v = ±h~ 1 (v v ± 1 -v u ), 



h , 

2^ 
h(D + 



f£>-), 

L>_) = h 2 D + D_ 



(5) 

Note that the operators So and p are dimensionless. 

Then, a direct but lengthy calculation starting from the definitions (|60|) leads 
us to the following expressions for D<>>") = D ( - 1 ' n '°'°'> , D^ n ~> = D& n >°<°\ the 
rest = (D^) 2 - and D^ 1 ' 71 ' 3 '^: 



n-l 



£)(!,«) 
R (n) 



h- l 8 Q + 



CkP 



k=l 



h~ 2 P 1 + J2 d k 



kP 



k=l 



n-l 



k=0 
s-1 



n + 1 



s-l 



k=l 



k=l 



where 



ak 



c.k = (-1 

s-1 



{2k + iy: 

{-iyc 2k 



kP p 



(6) 



j+k 



dk = kTT> 
= (-i) n+s 



r<2k+\ 

U s+fc 



2(n+l + fc)C 2 "+ s 



(7) 



Notice that the coefficients Ck and c?fe do not depend on n, while ak, bk do 
depend on n and s. 

The leading order truncation error of order 2n is defined as 



dv 
d.v 



x 

cPv 
d 2 x 



j-^(2,n) _ rp(2,n 



d 2n+1 v 



d 2n+2 V 



,2n 



rfa ,2n+2 



+ 0(/l 2,i+i ). 



h 2 " + 0(/i 2 "+ 2 ) . 
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| j^(l,n 7 s) | 


p 
i — 


n=2 


n=3 


n=4 


s=0 


i 

6 


1 

30 


1 

140 


l 

630 


s=l 


1 
3 


1 

20 


1 

105 


1 

504 



Table 1: Leading order truncation errors for the first order discrete derivative. 



A direct calculation yields 

ji(l,n,«,l) _ m(l,n,s,-l) _ ji(l,n,s) I i y+n ( n + s )K n ~ s )' 

1 ' (2n+l)! 

2(n!) 2 
(2n + 2)! ' 



= d„ = (-l) 



It is well known that the centered FDO has the smallest leading order truncation 
error, |T( 1,n '°)| < \T^^\ for s > (see also table [J). 

3.2 Fourier representation of difference operators 

We assume a finite grid defined by a set of N points, 

SJN) = {x v = uh, with h = 2n/N, Mv = 0, . . . , N - 1} , (8) 
and consider periodic grid functions v v = v mo( ir v m, decomposed as 

v v = v(w)b v {iS), (9) 

^eSoA,N) 

where 

K(lu) = (2n)- 1 l 2 e^ x " , (10) 

f {-7V/2 + l,...,iV/2}, if N is even 
°^ ; \ {-(iV-l)/2,...,(iV-l)/2}, if N is odd. { ' 

The set Su_{N) represents the set of discrete wave numbers, and in the space of 
periodic grid functions the set {b v (u>), uj £ S U (N)} forms a orthonormal basis 
with respect to the scalar product and the associated norm 

(v,u) h = ^2 ^ (x v )u(x v )Vh, \\v\\ h = (v,v) h . (12) 

with Vh = h. The quantities v{uS) represent the discrete Fourier coefficients. 
The scalar product satisfies the Parseval relation: 

(v,u) h = ^2 v*(u))u(u). (13) 
we5„(iV) 
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Let £ = ujhe Se(N) with 

f {-ir + 2n/N,...,ir}, if N is even 
i { ' \ {-ir + n/N, ...,tt->k/N}, if N is odd " 1 > 

Now apply the shift operator S k on a basis vector b v (w). This leads to 

S k e iux " = S k {£,)e wx % with S k {0 = e^ k . 

The function represents the discrete Fourier symbol of the shift operator. 

For any discrete operator D = J2k ak {h)S k the Fourier symbol is defined by 

De iu,x„ = D{^h)e iux ", with D(£;h) = ^ a k (h)S k (0 , 

k 

and for a general finite difference operator £)( m >™> s > e ) the symbol is 

n+es 

f){m, n ,s,e)^. h)=h - m J- j mmk e* k . (15) 

k— — n-\-cs 

For the elementary discrete operators (O we obtain 

D±($;h) = ±h- 1 (e ±i t-l), 

MO = where (5(0 = sin£, 

P(0 = -O 2 (0, where 0(0 = 2 sin |, 

and it is useful to note that £>+(£; /i) = £>_(£; ^) = /i _1 f2(0- 

The symbols for the first and second order derivative operators are straight- 
forwardly computed using ([6]), 

£>V>rifch) = -al^ n \i)h- 2 , 
RMfch) = r (n) (£)fr- 2 
D^'^fch) = (ed^' n ' s \(,) +id ( - 1 ' n ^{0) h- 1 , (16) 
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where we define 

n-l 



2 k- 



fe=0 
n-l 



f (n) = _(^«))a + ^) = lW ^1) ^ J^L^k j 

fe=0 n + 1 + ^' 

s-1 

j(i,„, s) = f2 2 "+ 2 ^(-l)"+ fe 6 fe n 2fe , 

fe=0 

s-1 

£l,n,s) = £l,n) + SCl 2n Y^(-l) n+k a k Cl 2k . (17) 

k=0 

In the following we list a series of particularly relevant properties of the Fourier 
symbols, further properties are given in appendix [Cl 

First note that the quantities f n and cft 2 ^ are positive, and even more, from 
p7|) it is easy to check that the following inequalities hold: 

< d (2 ' n) < S 2 < n+ V , (18) 

i>~: — >^ . (19) 

~ d 2 ^ ~ d( 2 .™+!) 

n-l 

C*,; 1 ^ 2 < il 2 < d^ n) < C n Cl 2 , where C„ = 1 + ^ \d k \ A k > 1 .(20) 

fe=i 

The real part of the Fourier symbol of the hrst derivative is an even function of 
the frequency £, while the imaginary part is an odd function. The real part of 
the Fourier symbol of the first derivative also 

• vanishes for centered operators (s = 0), 

• keeps the same sign for all frequencies, in case the operator is one-point 
off-centered (s = 1), 

• changes sign for off-centering by more than one point (s > 1). 
The derivatives with respect to £ of the Fourier functions satisfy: 

d^ 2 ^ = 2& 1 ' n) , <%f(™' = 2-^5-fi 2 "d( 1 ' n ) (21) 



(2n) 



^d (1 ' n,s) - ( -^r sin( S , 'V"' ~ = 1 - 7^ cos( S ■ 

C 2n L/ 2n 

In the following we show some plots to illustrate how the errors of the Fourier 
symbols scale with the order of approximation and off-centering. The error is 
defined in respect to the continuum limit, i.e., /i~ m (i£) m for i)( m >™> s > £ ). 



Figure [TJ shows the Fourier symbols S- X ' n \ <i( 2 <™) and f^ n > as functions of the 
frequency £ for different orders of accuracy. For increasing approximation order, 
the second derivative becomes more accurate for all frequencies, while the first 
derivative does not converge to the continuum limit for the highest frequency 
in the grid, where the symbol is zero. The 7r-frequency will not be captured 
also by the off-centered discrete operators associated with the first derivative. 
In addition, for them, the error scales with the order only at small frequencies. 

Figure [5] shows the scaling of the error for S- l ' n,s > with the off-centering at 
fixed order of approximation. In the region of small frequencies, off-centering 
increases the error. At larger frequencies this behavior changes. For each s, 
there are exactly s frequencies in (0, 7r) where the error cancels. However, for 
s > 2 there are large intervals where the error overcomes by far the error when 
s — 0. For s = 1 we observe that while at small frequencies, the error is slightly 
larger than for s = 0, for each order 2n, there is a frequency, beyond which 
the error is smaller than for the case s = 0. This frequency can be computed 
numerically, e.g. = 1.3787, £ (2) = 1.0036, £ (3) = 0.8234, £<4> = 0.7136. 




Figure 1: Fourier symbols as functions of the frequency £, for different ap- 
proximation orders. For increasing order the second derivative becomes more 
accurate for all frequencies, while the first derivative does not converge for £ = tt. 



3.3 Generalization to rf-dimensions 

In this work we will use the common straightforward generalization of finite 
difference operators from one to d dimensions. We extend first derivatives in a 
particular coordinate direction in the trivial way, and second derivative opera- 
tors in the ^-directions are defined as 



„<*»>_ J Df^D^\ for jtl 
°« - { Df'"\ for j = l. (22) 
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Figure 2: Errors for d^ 1,n,s ' at fixed order but different off-centerings arc 

shown, scaled with |c„|£ 2 ™ (n = 1,2,3,4). For small £ the curves are straight 
lines with the slope \c n T( 1 ' n > s > | , and larger s increases the error. At higher 
frequencies this behavior changes, but for s > 2 large intervals appear where 
the error overcomes by far the error when s — 0. For s = 1 we observe that 
while at small frequencies, the error is slightly larger than for s = 0, for each 
order there is a frequency beyond which the error is smaller than for s = 0. 

The Fourier symbols for the first and second derivative operators take the form 

bf^ = / ti 1 ( ( ,d( 1 ^)fe) + v:<i (, ^ ) fc)): 
bf' n) = hj\df' n \ 

[ ; ^ i 

Df n) = hjXU , (23) 

1 -dj* j=l 

In order to simplify notation, we will also use the following convention: any 
function in frequencies, and possible, grid spacings, /(^ , . . . & r ; , . . . hi r ) will 
be referred to by f^,,,^, and by / in case it depends on all frequencies. More 
detailed definitions for the d-dimensional case are given in appendix [B] 

3.4 Dissipation Operators 

In order to achieve numerical stability for problems that go beyond the linear 
constant coefficient case, it is common practice to add artificial dissipation to 
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the right-hand-sides of the time evolution equations. In this work we only deal 
with the constant coefficient problem, but in 15.21 we will also use dissipation to 
stabilize numerical schemes which would be unstable otherwise. 

Dissipation terms are typically chosen to converge away fast enough so as not 
to change the convergence order of the scheme. Here we use the Kreiss-Oligcr 
dissipation operator X>( 2m ) of order 2m [T] and its Fourier representation 2?( 2m ), 



Ti(2m) ( I)" 1 \ ^ „ i2ra-l/ n \«Vr> \ m -ft (2m) \ \ "* °J n2m 



d 

(24) 

for a 2m — 2 accurate scheme, where the parameters cr, > regulate the strength 
of the dissipation. Using this form of numerical dissipation, theorems can be 
proved concerning the numerical stability of non-constant-cocfficient hyperbolic 
PDEs pp. Note that it is more common to have the dissipation parameters <jj 
not depend on the direction or other parameters of the system. 

4 Numerical stability for first order in time, sec- 
ond order in space hyperbolic systems 

We now turn to the analysis of numerical stability for the system (J3]) , following 
[2] . This problem is greatly simplified by adopting the mcthod-of- lines approach 
where initially time is kept continuous and only space is discrctized (i.e. the 
semi-discrete problem). Then the discrete system to be analyzed becomes 

iv = Pv, v = {U 1 V) T 1 
at 



P = ( AW^+B C \ 

\ D 3l Df n) + EWf' n) + F (PDf' n) + J J { ' 

We consider periodic grid functions in each direction, and Fourier transform the 
system as discussed in appendix [EJ Then a first order reduction is performed 
by introducing the variable w, 

d 2 

w = ifl u, n 2 = \D+j , (26) 
j'=i 

where D + j is the Fourier symbol of the forward finite difference operator in 
the j-direction, D + j, The case fio = (which corresponds to zero frequencies 
in all directions) does not play any role in the stability analysis^, so we define 
5|(iV) = S ( {N) - d and assume £ G S^(N). 



2 the zero frequency vector corresponds to a term constant in space. 
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By (f2l))) we obtain the following system of ODEs 



— vr = PrVr with Vr = (u,w,v) 
at 



T 



Pe 



I B (ifl )-^Awf' n) C \ 

Aibf' n) + b m Q c 

\-i ( nil n( 2 '") _l T?j ni n( 1 >'0 



I F (iflo)- 1 {p»D)i' n > + EW) x ' nj J G J D y j ' ' + Jj 



.(27) 



Using the theorem 5.1.2 of [T] CHH show that the terms which correspond to 
the continuum lower order terms can be dropped from Pr without affecting the 
stability analysis if 

(iQo)- 1 ^, kbf> n \ kn^bf^ (28) 

are bounded for all frequencies £ £ S£(N). In the relations (|2"5]l . k represents 
the time step. We will show in lemma |4~T1 that this is indeed the case for any 
order of accuracy 2n. 

Having proved this, the rest of the discussion in CHH applies. The problem 
now reduces to the analysis of a first order system with the principal part: 



{itt^D^bf^ &bf' n) 



(29) 



For this type of system, sufficient conditions for stability have been deduced 
in pp. These conditions have been exploited in CHH to analyze the stability 
of the second order system. By introducing the so-called second-order principal 
symbol of the semi-discrete system, 

( AW [1 - n) C \ 
P ' = y Dj l£(2,n) Gjf) (l,n) j - (30) 

and assuming that the time integration is done using one-step explicit schemes, 
CHH show that the following conditions are sufficient for stability: 

Condition 1: There exists a hermitian matrix H (£, h) such that 

#~% < H < KI Qo , I Qo = diag[fig, I qN ], 

HP' + P n H = 0, (31) 

for some positive constant K. 

Condition 2: The eigenvalues of kP' have non-positive real parts and 

cr(fcP') < a (32) 

where a(kP') is the maximum spectral radius of kP' and ao is a constant specific 
to the time integrator. 
Remarks: 
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• The condition (|3"Tj) implies that the semi-discrete problem is stable with 
respect to the norms D± defined as: 



<^ ± =E W D * U t + \\V\\ 2 h . (33) 



where ||.||^ is the d-dimensional analog of (fT2|) . 

• The semi-discrete problem is stable also in the norm ||u||^ H defined by: 

\\< >H = E ( 34 ) 

This norm is conserved by the principal symbol of the evolution system, 
\Wt,.)\\ h>H = \\v{Q,.)\\ htH . 

• The constant ao in (|32[) denotes the radius of local stability on the imagi- 
nary axis (Ri s ia) in case the eigenvalues of kP' are purely imaginary, and 
the radius of local stability {Ri s ) 1 otherwise. 

• In case all the grid spacings arc equal, h\ = ■ ■ ■ = = h, and we introduce 
the Courant factor A = k/h, then the relation (|32[) provides the Courant 
limit: 

A<^_. (35) 
" a(hP') 

• If the right hand side of the system ([25| is modified by adding artificial 
dissipation (using the operator 2?( 2m ) defined in ([24| ) and/or by adding 
shift advection terms of the form Ifi^D^ 71 ' 83 (where ]jy' n,s " e ^ [ s the 
non-centered FDO in the j-dircction constructed from ©), these modifi- 
cations only have effect on the diagonal entries of the principal part. The 
new system will have different eigenvalues than P' but the same set of 
eigenvectors. The symmctrizcr will not depend on the way we discrctizc 
the advection terms, nor on the dissipation operator. The stability Con- 
ditions 1-2 remain valid if 

(ffio)- 1 ^ 1 '"'^, (zfto)- 1 ^ 2 ") (36) 

are bounded and this will be shown below together with the boundedness 
of the terms 



Lemma 4.1 The following quantities are bounded for all frequencies £ G S^{N). 
(iOc)" 1 ^, kfjf' n) , kSl^Df^, (iOo)- 1 ^ 1 '"' 8 ^, ^o)- 1 ^ 2 " 1 ) (37) 



3 for the classical fourth order Rungc-Kutta, Ri 3 i a = vo = 2.83 and Ri 3 = 2.61. 
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Making use of the relations (j2"3"]l and ([24"]). the proof reduces to showing the 
boundedness of 

From the relations (fl~7)) we observe that each of these quantities can be written 
formally as a product 1 fljF(^j), with F(^j) a continuous and bounded func- 
tion in (— 7r,7r]. Since f^fl,- is bounded for all Clj g (—2,2], j = l,...,d but 
not all zero in the same time, we obtain the desired result. 



5 Application: Scalar Wave Equation 
5.1 Semi-discrete Problem 

The system © is discretized assuming, for simplicity, that the grid spacings are 
equal (hi = ■ ■ ■ = h ( i = h). The case hi ^ hj for some directions i and j does 
not introduce further complications in the following analysis. 

We construct the semi-discrete system corresponding to (|2|) by: 

, d 
dt ^ 3 

3 = 1 

j d 

jK = -'/;;" <!> • Y > J ^! "" " (38) 

This way of discretizing the first order derivative terms, which correspond to 
advection along the shift vector (3 l , with off-centered derivatives has become 
customary in numerical relativity (see e.g. [Tl l27l[2"S] ). 
We define the shorthand quantity A as 



\ 



il j(l> n ) jC 1 )™) i ii-i n ) 

^ l dj 'dj + / J l n r) 

3 = 1 



Then the discrete symbol, the diagonalizing matrix and the eigenvalues can be 
written as 

/ «j£)(bn,s,e) l \ 

f ~ l= ( % ! ) ' A± = /^J 1 '"'"'"' ± ^ ■ WO) 

Because fj > 0, according to (JT7J) , and the matrix 7 j7 is positive definite, the 
quantity A is real and A > with equality only when all £j are zero. Thus 
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is a symmetrizer for the system ([55)1 . We observe that the symmetrizer does 
not depend on the diagonal entries of the symbol P' , e.g. does not depend on 
the way we advect the shift terms. 

We still have to prove that there exists a constant K > 1 such that 

K~ x £ll < A 2 < KUl . (41) 

The positivity of the matrix 7 j/ implies the existence of a constant c\ > such 
that 

ci < min t" and 7 J V j2 /, > ci \y | 2 , Vy = (y u . . . , y d ) € R d . (42) 
Furthermore, because 7 J '' < oo there also exists a constant c 2 > such that 

c 2 > maxf 3 and 7^3'^ < c 2 \y\ 2 , = (j/i, . . . , y d ) £ R d . (43) 
Using (|4"2"j) and the inequalities (f2"U)) we obtain 

h 2 A 2 > (mm^)J2 ?f ] + ci St 4 ? 1 '" 5 ) 2 ^ Cl E 4 2 '" } ^ Cl ^o ■ 

3=1 3=1 3=1 

On the other hand, by (|4"3")l and again (f2"0)l we have that 

/i 2 A 2 < (max 7 «) ^ fj n) + c 2 J2(df> n) ) 2 < c 2 ]T df n) < c 2 C n Cl 2 ■ 

3=1 3=1 3=1 

We chose K = maxjc^ 1 , (c2C„), 1} and obtain the relation (|4"Tj) . 

The conserved discrete quantity in physical space associated to H, i.e. the 
norm \\v\\ h H defined in (f34|) . is 



\ v fh,H = 



3=1 fc=l 37^ 



where v = (§ T , K T ) T . Having proved the existence of a symmetrizer we have 
proved that the semi-discrete problem is stable with respect to the norms D + 
and H . Note again that the stability property does in particular not depend on 
how the shift terms arc discrctizcd. 



5.2 Courant Limits and the Role of Dissipation 

In order for the fully discrete problem to be stable we impose the non-positivity 
condition on the real part of the eigenvalues and restrict the Courant factor A 
according to the inequality (1351) : 

i?e(A±) < 0, (44) 







max 


hA± 


tes L 





(45) 



15 



From (|4"Uj) we have 

d 

W?e(A±) = '"'^ . ( 46 ) 

3=1 

The relation (|4"4")) has to hold for all frequencies £j E (—tt,tt). Because each 
term j in the sum (|46p can be canceled individually at £j = 0, the non-positivity 
condition has to applied for each term. The problem reduces to the study of 
the one-dimensional case, 

/3ed (1 ' n ' s) (£) < with £e(-7r,7r] (47) 

Because dl 1 '"'*' is zero for s = 0, negative for s = 1 and changes sign for 
s > 2, it is clear that the condition holds for centered and one-point upwindcd 
(e = sign (3) schemes and is violated in all the other cases. 

However, the condition can be reestablished if appropriate artificial dissipa- 
tion is added to the system. By using the Kreiss-Oligcr dissipation operator 
(j2"4")h the condition (|4T|) changes to 

^d(i^)(^)__i_ (7 f 2 2(„ + i) (o < with ^(.^ (48) 

This imposes a lower limit on the dissipation parameter a: 

{ 2 2(«+i) \p\ & M ) e = sign p ( U p W i n d ) 
(49) 
2 2(n+i) ^1 g(_">s) ) e = _ sign p (downwind) 

where we have denoted 

, \ H(l,™>s) i \ AC 1 .".*) 

a! = max — , o_ — — mm — . 51) 

+ f2e(0,2]fi 2 («+ 1 ) oe(o,2]ft 2 (™+ 1 ) 

and used the fact that d( 1,ra ' s ) is, according to (JTTJ) , a sum over powers of Cl. 
In table [2] we give the formulas for cr^™ for s = 1, 2, 3. We remark that for 

all 7i > 1, < 0, a { "' 1] > and a ( ±' s) > for all s > 1. 

This means that when using one-point upwinded stencils, we can add "nega- 
tive" dissipation and still obtain a stable scheme. In fact, the following situations 
are equivalent: 

• Upwind one point and add dissipation with a = 2 2 (" +1 ) \/3\ Oj™' 1 ' < 0. 

• Downwind one point and add dissipation with a = 2 2 ( Tl+1 ) \p\ g-*™' 1 ) > o. 

• Use the CFDO operator 1/2 (zKb'W) + £)(i,n,s (constructed with 
2(n + s) + 1 points), and do not add dissipation, a = 0. 
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s=l 


l l 


1 1 


2C^+ 1 "+ 1 


2CJ+ 1 n+1 


s=2 


1 2 


1 2n 


2C n+2 n+1 


2C 2 "+ 2 n 2 +3n+2 


s=3 


1 ra(n+4) 


1 3 


2C r. + a ( n+ l)( n+ 2)-J 


2CJ+ 3 n+1 



Table 2: Formulas for the dissipation parameters a±' s ' when s = 1,2,3. The 
quantity 2 2 ( n+1 ) J/?- 7 1 a± (± stands for upwind/ downwind) represents the min- 
imum dissipation that one has to add to make the numerical scheme stable. 



In any of the above three situations, the real part of the eigenvalues is zero, so if 
the Courant limit is small enough then we obtain stability on the imaginary axis. 
Now, coming back to the d-dimensional case, it is obvious that the dissipation 
parameters in ([M]) . aj, have to be chosen according to the value of the shift 
and the type of off-centering in the j-direction (aj > cr m i a (/3j,n, Sj, tj)). It can 
shown that, by choosing exactly Oj = a m - ln (/3j ■, n, s 3 ■, €j) the Courant limit is 
maximized. However, to compute it explicitly, (as a function of shifts, order of 
approximation, and off-centerings) is not easy in the general case. 

In the particular case of a flat d+ 1-metric with zero shift, the Courant limit 
is easy to write down: 

- 2v^o:' 

where C n is given in (|20p and ao = Ri S i a - In the general case, the Courant limit 
has to be evaluated numerically. 

For the TD wave equation with shift (3 > with upwind discretization of 
the advection term and adding the minimal amount of dissipation if necessary, 
the limit of the Courant factor is given by 



A ( "' s) (/3) 



max 

£e(7T,7r] 



We compare the Courant limits for different orders of approximations at fixed 
advection stencil in Figure [31 and the Courant limit at fixed order of approxi- 
mation for different advection stencils in Figure |4j 

In Figure [3] we see that if s = 0, the higher the order of approximation, 
the lower the Courant limit. For s > 1, this is not true anymore beyond a 
certain value of the shift. For large shifts we observe that increasing the order 
of approximation, actually decreases the Courant limit. 

Comparing different stencils in Figure 21 we observe that advecting more 
points decreases the Courant limit, and there is a significant drop in the Courant 
factor between s = 1 and s = 2, for all orders of approximation. 
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Figure 3: Courant limit as a function of /3, for different orders of ap- 
proximation at fixed off-centering s. For s = (left) no dissipation is 
needed, (a = 0), and we are in the regime of local stability on the imaginary 
axis (ao = 2.83). For s — 1 (middle), again no dissipation is needed (<r = 0), 
but now we are in the regime of local stability (ao = 2.61). For s = 2 (right) 
dissipation is required and we add the minimum amount in order to attain 
stability. 




Figure 4: Courant limit as a function of /3for different advection stencils 
at fixed order of spatial accuracy. From left to right: Courant limits at 
approximation orders 2, 4, 6. As in Figure [3] the Courant limit calculation takes 
into account whether we are in the regime of local stability on the imaginary 
axis (the case s = 0), or only local stability (for s > 1), and the minimal amount 
of Kreiss-Oliger dissipation is added for s > 2. 
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5.3 Phase and Group Speeds 

For the wave equation the continuum phase and group speeds are: 

A , d 



v„ = v„ = B n ± v/7™, where v„ = — , v„ = n J — — A . 

ujq dujj 

where /3™ = ftrij and 7™" = j^rijUi. The discrete speeds can be defined in a 
simillar manner from the discrete eigenvalues A. This would lead to complex 
speeds in case of off-centered schemes which are difficult to investigate [25] . 
In the following, we assume that the real part of the discrete eigenvalues has 
been cancelled by adding appropriate artificial dissipation terms (positive or 
negative). Otherwise, our investigation remains valid only at small frequencies, 
where the damping/amplification effect introduced by the real part is dominated 
by the dispersion effect associated to the imaginary part of the eigenvalues. With 
these remarks, we define the discrete speeds by: 




± — Uj 



We also restrict attention to the one dimensional case. Because ± speeds inter- 
change when £ changes sign, it is enough to consider only the "+" speed over the 
whole spectrum £ e (— ir,ir}. Also because we will compare speeds at different 
orders of approximation or at different stencils, we attach the superscript (n, s) 
(or only (n) in case s = 0), to the symbols representing the discrete speeds and 
the corresponding errors: 

z)("' s) (0 = I (p$i.»,>) + v 7 ^ 2 '™)) , 

(£) = A + Vd(2,")) . 

The continuum limits for both, phase and group speeds are (3 + 1 for £ > and 
/3 — 1 for £ < 0. We will analyze the behavior of the speed errors defined as 

- + (51) 

4"'" = fi(^<i <1 '"''>-l) + (^\ / 5M-.iign«) . (52) 
We will also assume j3 > without restricting generality, if /3 — > — /?, then 
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5.3.1 Small Frequencies 

When £ ~ one can show that the phase and group speed errors satisfy 



(n!) 2 r 2(n + l) 

4"< s > = -(2n + l)|c„ f 



( 1 ^ (rt + s)!(n-s)! ^ g , sign £ 



(™!) 2 ^ 2(n + l) 



e 2,l + o(^ 2 ' i+2 ), 

e 2n + o(e 2n+2 ) 



Because the errors scale with £ 2n , it is obvious that for small enough frequencies 
higher order approximations will improve the phase and group errors for all the 
values of the shift and for all advection stencils. 

If we keep the order fixed and compare the speeds corresponding to an off- 
centering by s > 1-points with the ones corresponding to the centered scheme, 
s = 0, then one can easily show that the off-centered scheme improves over the 
centered one 

• the "+" numerical speeds (£ > 0) if s is odd and (3 is small enough 

• the "-" numerical speeds (£ < 0) if s is even and j3 is small enough 
where small enough means 

1 1 

/?< 



( n + l) (n+ s y.(n-s)\ _ , ■ 



Obs. For s = 1, the above inequality becomes j3 < > " . Also notice that with 
increasing s the above limit on j3 decreases. 

In the next subsection we will analyze the behavior for the whole spectrum in 
some more detail. 

5.3.2 Comparison with wave equation written in first order form 

If the wave equation is written in first order form (approximating the hrst deriva- 
tives with the corresponding CFDO), then the eigenvalues become (/iA±)(£) = 
i (/3 ± 1) dP' n \ For £ ~ one then gets 

4 n) = -03+ sign o\c n \e n +o(e n+2 ) 

4 K) = - W + «gn 0(2n+l)\c n \e n + O(e n+2 )- 

We notice that for a given order, the second order system discretized with 
CFDO, has smaller phase and group errors then the first order one (for both 
eigenvalues), if and only if \j3\ < gT^+fy- If 1/3 1 is not in this interval then one 
pair of speeds (phase and group) is better approximated by the second order 
system, while the other one is better approximated by the first order system. 
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5.3.3 Scaling of the Speed Errors with the Order of Approximation 

Lemma 5.1 If/3 = 0, then increasing the order of approximation decreases the 
phase and group speed errors for all frequencies. 



To prove this we make use of the relations (|2lj) in the definitions of the speeds 
and obtain 



-,(«) 



#l,n) 



;(n) 



sign £ , 



- sign £ . 



Using the inequalities (fT5]l and p^|) one can easily show that 



i(n+i) 



< 



-(«) 



and 



< 



for all frequencies. The situation is illustrated in Figured] 



where we plot the speeds and Vg' l> versus £. 



Phase Speed when j6=0 
1 



Group Speed when /J=0 
1 




Figure 5: Phase and Group Speeds for /3 = 0. The higher the order of 
the approximation, the more accurate the phase and group speeds are for all 
frequencies. 

If /3 7^ then it is not true anymore that higher order approximations im- 
prove the numerical speeds for all frequencies (not even for the case of using 
CFDO). Though one can go into details and determine the regions in the spec- 
trum where the scaling with order fails, we restrict ourselves to illustrating 
this situation by plotting the numerical speeds versus frequency at a particular 
value of the shift. In Figure [6] we show the numerical speeds at different orders 
of approximation with the same advection stencil when /3 = 0.5. 



5.3.4 Scaling of the Speed Errors with Off-centering 

The next question we want to answer is what happens with the numerical speed 
errors, if we keep the order of approximation fixed and vary the off-centering of 
the first derivative. For example, we illustrate this situation in Figure [3 when 
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Figure 6: The phase and group speeds at different orders of approximation but 
keeping the level of off-centering fixed, for /? = 0.5. 

/3 = 0.5. As one can see already in these plots, although we know that off- 
centering increases the error of the finite difference operator, it is not necessary 
that the numerical speeds will follow the same pattern. E.g. for this value of 
the shift, the "+" speed seems more accurate with s = 1 than with s = for 
all the spectrum. In the following we will determine the regions in the (£, /3)- 
plane where off-centering improves the numerical speed errors over the centered 
scheme. 



Phase Speeds 



Imposing 



where 



4n,s) 



< 



Jn,0) 



and using the definition (|51[) yields the inequality 



(53) 



f (n,s ){0 



( J(l,n 1 «) + ^(l,n I 0)_2^ j 

5 (n) (0 



(54) 



We have g( n \C) > but can change sign over the spectrum. The inequal- 



ity ([531) holds at a given frequency £, if /3 > ^™> s )(£) and sign fffliO < 
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Figure 7: The phase and group speeds errors (scaled with £ 2n ) are shown at 
different advection stencils for f3 = 0.5. 



or /3 < (3( n ' s \£) and sign (£) > 0. In general, the regions in (£,/?) plane 
where at fixed order of approximation, off-centering by s points improves the ac- 
curacy of the phase speed, are difficult to determine analytically and we restrict 
ourselves to a numerical evaluation (Figure [8|) . What we see in the plots is that 
if s is odd (even) then for sufficiently small /3, the "+" ("-") speed has smaller 
error compared with the case of CFDO in some intervals of the spectrum that 
include the small frequency range. However these regions become narrower with 
increasing the off-centering, such that for s = 1 we have the strongest effect. 
We analyze this case in more detail below. 
If s = 1, the functions fi2^ and /3^- n,s \ defined in |54|) become 



/i (n,1) (0 



\Cn-l\ 



(sin£)n 



2 It 



(55) 



We have sign /i (£) = sign £ and fi' l ' 1> (±n) = 
Then the inequality e p ™' 1 '' < 



(n,0) 



holds 



• for £ > if ^"'^(O < or < (3 < /3 (n ' 1) (0> 

• for £<0 if /3 > P (n > 1} {0 > 0. 



23 



The limits of /3 (ti,x) (£) in and tt are 



= -lim/^"' 1 ^) 



1 ' 



-1 + 2- 



<0. 



It can be shown that the equation (3 = /S^ 1 ' 1 '^) has at most one solution in 
each of the branches £ > and £ < 0, that we will denote by £ ± . It turns out 



that 



-(n,i) 



< 



=(n,0) 



holds if 



• p < 1 

• 1 -2^ 



^ and £e (o,7r), 

</3< ^ and£e(-7r,r)U(0,7r), 



• /3>^T and£G (-7r,r)U(£+,7r). 

At a given order of approximation, 2?t,, for sufficiently small f3 < -^j-^ the "+" 
speed has smaller error in the case when we advect one point than in the case 
when we use CFDO, for all frequencies < £ < 7r, but the "-" speed will have 
larger error, at least for small and mid frequencies. 

If /3 > then for both ± speeds, in the regime of small frequencies, 

the CFDO give less error than one-point advected scheme, while for mid and 
high frequencies the situation reverses. The interval of small frequencies where 
CFDO are better than advected scheme shrinks with increasing the order of 
approximation. 



Group Speeds 



Imposing 



where 



< 



-(n,0) 



and using the definition (|52[) yields the inequality 

£(«.»>(£)) 



Fi („, s)(0F K s ) (0 ^ 



< 0. 



p(n,s 


(0 






p (ra,s 
r 2 


(0 






G (n 


(0 


= d^ n 






(0 


G(") 

b 2 


(0 

'(C)' 



(56) 



and/{" 2 ' sJ and ff (") are given by d54j) . It is easy to see that G (n) (£) = -G (n) (-£)- 

However the signs of F^ 8 \^) are more difficult to determine. As in the case 
of phase speeds analysis, we determine graphically (see Figure [9]) the regions 
in (£, /?) plane where at fixed order of approximation, off-centering by s points 
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Figure 8: Shown are the regions where advectcd stencils improve the phase 
speed error over the centered scheme. The regions are delimited by the quantity 
/3(™> s ) and the zeros of the function f( n ' s > as defined in ([5l 



improves the accuracy of the group speed. We see the same qualitative behavior 
as for the phase speeds, in the sense that for sufficiently small /3, the "+" ("- 
") speed has smaller error compared with the case of CFDO at least at small 
frequencies, and off-centering decreases the extent of these regions in (£, 
space. 
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In case s = 1, the relations ([55)1 become 



^(0 = (n + 1 ^ |c "- l| (^ I + cosg) n»" 

^)(0 = (n+1 ^ |C "- ll (-^+cose)^ 



2 (sign £ - dC 1 .") / Vd( 2 ' n ) ) 
/S (n,1) (0 = — -• (57) 



By analyzing the monotony of these functions using the properties from l3.2I the 
following result can be formulated: 

At a given order of approximation 2n, for sufficiently small j3 < ^h;, the 
"+" group speed has smaller error in the case when we advect one point than 
in the case when we use CFDO for all frequencies < £ < n — arccos ^py, (in 
the case of phase speed this was the whole range (0,7r)!), but the "-" speed will 
have larger error, at least for small and mid frequencies. 

If /3 > ^pj- then for both ± speeds, in the regime of small frequencies, 
the CFDO give less error than one-point advected scheme, while for mid and 
high frequencies, the situation reverses. The interval of small frequencies where 
CFDO are better than advected scheme narrows with increasing the order of 
approximation. 



5.3.5 Centered versus One-Point Upwinded Scheme, Numerically 

In l5.3.4l we showed that when < j3 < the numerical "+" speeds are better 
approximated with one-point off-centered schemes than with centered schemes 
at least up to very high frequencies in the grid. In this section we show some 
simple numerical tests to illustrate this fact. We chose Z-periodic initial data: 

$(o,i) = e -( 2 ^ 2 r lsi " 2 (f*-f) 5 

K(0,x) = ad x $(0,x), x€[0,l). (58) 
The parameter a E [—1,1] sets the amplitude of the "±" components, 

C ± = {a±l)d x <t>. 

When a = 1(— 1) the signal is purely "left" ("right") going and when a = 0, the 
signal is equally distributed between both modes. 

We choose a grid with N = 101 points and resolution h — 0.01, the width 
of the grid is I = Nh = 1.01. Also we chose r = 0.1, j3 = 0.5 and we integrate 
the wave equation using fourth order FDOs for space derivatives and the fourth 
order Runge-Kutta as time integrator. 

We let a £ {1,0,— 1} and for each value of a we look at the errors for the 
main variables when s = 0, 1 (see Figure [TP)) . The numerical results show that, 
indeed, when the signal is "left" going, the upwinded scheme has less error than 
the centered scheme, while when the signal is going "right" , the centered scheme 
is to be preferred. 
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Figure 9: Shown are the regions where advected stencils improve the group 
speed error over the centered scheme. The regions are delimited by the quantity 
/3(™> s ) and the zeros of the function F[ n ' s ^ defined in (15F 



6 Conclusions 

In this paper we have investigated several aspects related to the discretization of 
the initial value problem for first order in time and second order in space systems 
of differential equations, using high order finite difference operators. Special 
attention has been paid to the situation when some of the first derivatives are 
approximated with off-centered discrete operators as is customary for treating 
black hole spacetimes in numerical relativity. Our investigation has been divided 
into three parts: (a) We started with an analysis of certain properties of the finite 
difference operators (Scction[3]). (b) Using these properties we have extended the 
validity of an existing stability method (Section^]), (c) We analyzed the stability 
and the numerical speeds in the case of the scalar wave equation (Section [5]). 
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Figure 10: The top panel shows the errors in the I 2 norm, the bottom panel 
snapshots of $ at t = 99CT. The three cases from left to right are a purely 
left- going signal (a = 1), a signal with equal amplitudes of left and right going 
modes (a = 0), and a purely right-going signal (a = —1). Red lines mark the 
centered scheme, green the one-point advected stencil, blue the exact solution. 
For a = 1 upwinding is more accurate, but the situation reverses if the signal is 
right-going, for a ~ both schemes yield similar accuracy. 
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In the following we will give a brief overview of the results. 

(a) Analysis of first and second order discrete derivative operators 

A set of mathematical properties have been deduced for the Fourier symbols 
associated with the second order centered and first order (not necessarily cen- 
tered) discrete derivatives. They are in the form of inequalities, recursive and 
differential relations for the Fourier symbols at different orders of approximation 
or at different off-centerings. Here we mention two of them: 

While first derivatives do not converge in the limit n — > oo at the maximum 
grid frequency (£ = tt), second derivatives do converge at all frequencies (that 
is the highest frequency in the grid will not be captured by the first order 
derivative, regardless of the order of approximation or the off-centering, while 
the second centered derivative can "see" it and approximates it better with 
increasing order). 

For first order derivatives, increasing the off-centering (s) at a fixed order 
of approximation increases the error of the derivative at small frequencies. At 
larger frequencies this behavior changes. E.g., for s = 1, beyond a certain 
frequency £' n \ the error is smaller than for the case s = 0. As a consequence, 
off-centering of the first order derivative in the case of the advection equation, 
increases the error at small frequencies, while at high frequencies, this situation 
can reverse. 

(b) Generalization of a stability analysis method 

In [3], necessary and sufficient conditions for stability have been deduced as- 
suming that ([3]) is discretized using 2nd or 4th order CFDOs and integrated in 
time using a time integrator locally stable on the imaginary axis. The validity 
of this stability method is extended here to 2n-order spatial accuracy, including 
also the case where some derivatives arc approximated with non-centered FDOs 
and dissipation is added to the system. It is pointed out that neither adding 
artificial dissipation nor shift advection terms affects the eigenvectors of the dis- 
crete symmetrizcr, and thus the conditions for semi-discrete numerical stability. 
The Courant limit will of course be affected in general. 

(c) Application: Scalar Shifted Wave Equation 

The stability method presented in Section |4] is applied to the case of the wave 
equation on a curved background in d + 1 dimensions. 

In the case of 1-D shifted wave equation in flat spacetime, the Courant limits 
and the numerical speeds have been analyzed in detail in respect to the order 
of approximation and off-centering of the first derivative. 

• Courant limits 

Off-centerings by more than one point require dissipation for stability. In these 
cases, the minimal Kreiss-Oliger dissipation needed for stability has been com- 
puted and found to be proportional to the shift j3. For centered schemes, higher 
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order approximations have lower Courant limits. Interestingly, this does not 
hold for off-centered schemes (when adding just dissipation to be in the lo- 
cal stability regime) — for large enough shift, the Courant limit is actually 
larger for higher order schemes. Off-centering generally reduces the Courant 
limit drastically, except for at least fourth order accurate schemes, when only 
one-point off-centering is used: for higher than fourth order schemes one-point 
off-centering only leads to a minor reduction of the CFL factor. 

• Numerical speeds 

Without shift, higher order approximations always result in more accurate nu- 
merical speeds, with nonzero shift this is not generally true at higher frequencies. 

Although the truncation error for the first order derivative increases with the 
off-centering, the mixing with the second order discrete derivative in the scheme, 
causes upwinded stencils to give a higher overall accuracy in some situations. 

More precisely, it is shown that advecting shift terms by an odd (even) 
number of points reduces the errors of the "+" ("— ") numerical speeds in some 
intervals of the spectrum that include the small frequency range, if the shift 
is not too large. The extent of the regions in the (frequency, shift)-paramctcr 
space where this improvement appears, decreases with off-centering, in such a 
way that for s = 1 one gets the strongest effect. 

Thus, at a given order In. if the shift satisfies < /3 < then off- 

centering by one point has in comparison with the centered scheme, better "+" 
phase speed error for all frequencies, and better "+" group speed error for all 
frequencies up to a very high frequency in the grid, 7r — arccos ^rj- 

If the wave equation is written in first order form and discretized using 
CFDO, then for a given order of approximation, the second order system dis- 
cretized also with CFDO has smaller phase and group errors than the first order 
one, if |/3| < jT^+fy- If |/?| is not in this interval then one pair of speeds (phase 
and group) is better approximated by the second order system, while the other 
one is better approximated by the first order system. 

A detailed understanding of finite difference algorithms for first order in 
time, second order in space systems, in particular as applied to the Einstein 
equations, will require significant further work. Already for the shifted wave 
equation, it will be interesting to study the errors in the multidimensional case, 
e.g. when the wave propagates in a direction that is not aligned with the grid. 
A similar analysis for the full Einstein equations will require a substantial use 
of computer algebra methods. We also point out that for the Einstein equations 
much of the complications come from the nonlinear source terms, which are 
beyond the scope of our present analysis. 
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A Explicit Expressions for Finite Difference Op- 
erators 

Explicit formulas for general finite difference operators in one dimension can be 
constructed in a surprisingly simple way by use of the in the following lemma. 

Lemma A.l Consider a general finite difference operator £K m >™> s ' e ) ; which is 
written as a linear combination of shift operators of the form O) 

n+es 

D (m,n.s,e) = ^-m £ / m , n ,,, e , fc S fc . 
fc=— n+es 

Then the coefficients f m ,n,s,e,k are the coefficients of y k in the Taylor expansion 
of the function 

f m ' n ' s ' e (y)=y n - es 0ny) m 

around the point yo = 1 up to the order (y — yo) 2n . In general, the accuracy of 
this operator will be 2n + 1 — m. 

To prove the lemma, it is enough to consider the scalar function v(x) = e luJX with 
w6C and the associated grid function. By applying accordingly the differential 
and discrete operators we get: 



D {m, n ,s,e) Vv = h -m ^ 



d m v{x v ) = {iuj) m e tu:hw 

n+es 

fm,n,s 7 e.k& 
k— — n-\-es 

d m \i(x v ) = D^ m ' n ' s '^v u + 0(h q ) (59) 
Introduce y = e which gives iui = h^ 1 \ny. The relations ([51?)) lead to 

y n -"{hiy) m = f™ mk y k+n - es + 0{hi+ m )y- v+n -™ . 

k— — n-\-es 

In the limit h — > 0, y — > 1. The function y n ~ es (\ny) m is now Taylor expanded 
around the point yo = 1 up to (y — yo) 2 " an d the coefficients f m ,n,s,e,k are 
identified. What remains is: 

0((y - l) 2 ^ 1 ) = 0{h q+m )y~ l ' +n -' LS . 

After replacing y with its definition and taking the limit h — > we obtain 
q = 2n + 1 — m and the lemma is proved. 
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Corollary A. 2 The FDOs associated with the first and second derivative are 



n+es 



D (l,n,s,e) = h -l J2 a ntStt>k S k , 
k— — n-\-es 

£)(l,n) _ £j(l,n,0,0) _ ^-ly Wnjk 

fe=l ^ 



D (2,n) _ £) (2,n,0,0) = ^ (S* + S~ k ) , (60) 

where 



fc=0 



(-l)* + 1 (n+«)!(n- fl )! , , n f 2(-l) fc + 1 (n!) 2 , , 

fc(ri+es-fc)!(n-es+fc)! > ^ / U k 2 (n+k)\(n-k)\ 

a«,s,e,fc = { and j3 n>k = < 

e(H n - s -H n+s ), k = I -ELi^fe fc = 0- 

In the relations above, H n = * s the harmonic number. Note that 

kf3 n , k = 2a„, ,o,fc f° r k>l. 



B Finite difference operators in ^-dimensions 

Consider a d-dimcnsional grid defined by the set of points x v = (vih\, . . . , Vdhd), 
where v — {v\, . . . , Vd) is a multiple index, i/j £ Z and hj represents the grid 
spacing in j-direction (j = l,...,d). Corresponding to the continuum vector- 
function v : M. d — >• C x • • ■ x C we associate the grid vector-function v such that 

The shift operator by fc-points in the j-direction, S k is defined by 

SjV{ Vl ,..., Vd ) = V( Vlt ... v . +k ,., Mi) . (61) 

A discrete operator Dj acting in the j-dircction is constructed as a linear 
combination of the shift operators defined in (|61[) using the same weights as 
the corresponding one dimensional operator D; an operator Dj lm . m j r acting in 
ji . . . j r -directions j s constructed as a composition of one-directional operators 

( D h D D : 

Dj = "^2a k {h 3 )S!j where D = ^ a k {h)S k (62) 

k k 

D h ... jr = <nh, ..../-, )!)) ...!)[ (63) 

In order to represent the functions in Fourier space, we consider only grid func- 
tion which are periodic in each direction and limit the grid to having a finite 
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number of points, Nj for the direction j, j = 1, . . . , d. We introduce 

,N d ), 



h 


= (hi,...,hd), 


N = (iVr, 


b v {w) 


= (2ir)- d / 2 e luJX 




v h 


= hi...h dl 




SJiN) 


= 5 £ (M) x • • • 


x SJiN d ), 


S„{N) 


= Su(Ni) x • • ■ 


x SJiNa), 




= S € (Ni)x--- 





In the relations above, Sx(Nj e N), SJ^Nj € N)) and <%(iVj e N) have been 
dchncd in ((HJ, (fTTj) and (fT4|) . respectively. Then, the formulas for the Fourier 
decomposition ((9]), scalar product and a norm (fT2|) and Parseval relation (|13|) 
are valid also in (^-dimensions. 

Let £ <G S$(N) and apply the shift operator by fc-points in the j-direction 

on a basis vector This leads to 

S*b v (cj) = § k (t j )b v (u), with S fe (£) = e l « fc . 

Then the Fourier symbols for the operators Dj and Dj 1 ...j T from (|62H63p are 
dchncd by: 

Djb v (u) = D(& h i) b A") 
D il ...i T .b v {uj) = £)(&!,... £i r ;h il ,...hi r )b v (cj) 
with 

D&ihj) = ^^(/li)^(ei), 

A; 

D(£ il ,...£ ir ;h il ,...h ir ) = o(/N 1> .../ii p )^ 1 tei)----D r (6r) 

C Further Properties of the Fourier symbols 

1. Recurrence relations: 

JU.n+l) = d {1 ' n) +5\c n \(l 2n . 
£2,n+l) = d(2,«) + | dn |02„+2_ 

#!,»,.) = Jd.n.a-l) (-1)' Q2n +1 cog (2S ~ IK 



(n + s)<T +s 



(n + S )C 2 "+ s 



2;; 
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2. Small frequency behavior: 



2n+ 1 



(-1)" 



1 - 



2n + 2 
(-1) 



rp(l,n,s)£2n+2 



3. Sums 



E \ck\x 

k=0 

oo 

k=0 



2k 



§ v 1 - (§) 



,2/,' 



4. Limits n — > oo: 



lim d (2 ' n) (0 
lim ^ X '"'^(0 

?l — >■ OG 

lim d^ n ' s ~>^) 

n— >oo 

lim f 

5. The D< 2 '"'-norm defined by 



Vfe[0,7r], 

e, vee[o,7r), 

o, vee[o j7 r], 

0, V£ 6 [0, tt) . 



^ d n 



i 2 



is equivalent with the D+ norm. This norm has been used to prove strong 
stability of the initial boundary value problem for the wave equation in |16j 
for the second and fourth order accuracy case. 
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